Load libraries and data
easypackages::libraries("here","ggplot2","caret","e1071","pheatmap","reshape2","NbClust","grid","patchwork","readxl","patchwork","WGCNA","psych","nlme")
source(here("code","ndar_functions.R"))
source(here("code","euaims_functions.R"))
source(here("code","get_ggColorHue.R"))
options(stringsAsFactors = FALSE)
# Z-score threshold to use for subtyping
z_thresh = 0.8
codepath = here("code")
datapath = here("data")
figpath = here("figures")
resultpath = here("results","ndar")
plotpath = here("plots","ndar")
# read in data
Dverbal_Discovery = read.csv(file.path(datapath,"tidy_verbal_disc.csv"))
Dverbal_Replication = read.csv(file.path(datapath,"tidy_verbal_rep.csv"))
make_subtype <- function(data2use, z_thresh, mean2use=NULL, sd2use=NULL){
# compute difference score
vars2use = c("dbaes_atotal","dbaes_btotal")
diff_score = data2use[,vars2use[1]] - data2use[,vars2use[2]]
# compute mean and sd if necessary
if (is.null(mean2use)){
mean2use = mean(diff_score)
} # if (is.null(mean2use))
if (is.null(sd2use)){
sd2use = sd(diff_score)
} # if (is.null(sd2use))
# compute z-score
data2use$z_ds = (diff_score - mean2use)/sd2use
# make subtype factor
data2use$z_ds_group = "SC_equal_RRB"
data2use$z_ds_group[data2use$z_ds>z_thresh] = "SC_over_RRB"
data2use$z_ds_group[data2use$z_ds<(z_thresh*-1)] = "RRB_over_SC"
data2use$z_ds_group = factor(data2use$z_ds_group)
return(data2use)
} # function make_subtype
vars2use = c("dbaes_atotal","dbaes_btotal")
# compute Discovery mean and sd
ds_disc = Dverbal_Discovery[,vars2use[1]] - Dverbal_Discovery[,vars2use[2]]
mean2use = mean(ds_disc)
sd2use = sd(ds_disc)
Dverbal_Discovery = make_subtype(data2use = Dverbal_Discovery,
z_thresh = z_thresh,
mean2use = mean2use,
sd2use = sd2use)
# compute Replication mean and sd
ds_rep = Dverbal_Replication[,vars2use[1]] - Dverbal_Replication[,vars2use[2]]
# mean2use = mean(ds_rep)
# sd2use = sd(ds_rep)
Dverbal_Replication = make_subtype(data2use = Dverbal_Replication,
z_thresh = z_thresh,
mean2use = mean2use,
sd2use = sd2use)
maxScores = c(3,4)
p_disc = ggplot(data = Dverbal_Discovery, aes(x = dbaes_atotal, y = dbaes_btotal, colour = z_ds_group)) + geom_point() + xlab("SC") + ylab("RRB") + ylim(0,1) + xlim(0,1) + ggtitle("NDAR Discovery")
p1_top_left = p_disc + guides(colour=FALSE)
ggsave(filename = file.path(plotpath, sprintf("final_NDAR_Disc_scatterplot_z%s.pdf",as.character(z_thresh))), plot = p1_top_left)
p_disc
table(Dverbal_Discovery$z_ds_group)
##
## RRB_over_SC SC_equal_RRB SC_over_RRB
## 178 538 173
cor_res = cor.test(Dverbal_Discovery$dbaes_atotal, Dverbal_Discovery$dbaes_btotal)
cor_res
##
## Pearson's product-moment correlation
##
## data: Dverbal_Discovery$dbaes_atotal and Dverbal_Discovery$dbaes_btotal
## t = 6.6829, df = 887, p-value = 4.134e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1554332 0.2806578
## sample estimates:
## cor
## 0.2189469
p_rep = ggplot(data = Dverbal_Replication, aes(x = dbaes_atotal, y = dbaes_btotal, colour = z_ds_group)) + geom_point() + xlab("SC") + ylab("RRB") + ylim(0,1) + xlim(0,1) + ggtitle("NDAR Replication")
p2_bottom_left = p_rep + guides(colour=FALSE)
ggsave(filename = file.path(plotpath, sprintf("final_NDAR_Rep_scatterplot_z%s.pdf",as.character(z_thresh))), plot = p2_bottom_left)
p_rep
table(Dverbal_Replication$z_ds_group)
##
## RRB_over_SC SC_equal_RRB SC_over_RRB
## 196 533 161
cor_res = cor.test(Dverbal_Replication$dbaes_atotal, Dverbal_Replication$dbaes_btotal)
cor_res
##
## Pearson's product-moment correlation
##
## data: Dverbal_Replication$dbaes_atotal and Dverbal_Replication$dbaes_btotal
## t = 7.1459, df = 888, p-value = 1.864e-12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1700824 0.2943934
## sample estimates:
## cor
## 0.2331903
# run validation
# make subtypes using z-scores computed from the mean and sd of the training set
train_data = Dverbal_Discovery
test_data = Dverbal_Replication
mean2use = mean(train_data[,vars2use[1]] - train_data[,vars2use[2]])
sd2use = sd(train_data[,vars2use[1]] - train_data[,vars2use[2]])
tmp_train = make_subtype(data2use = train_data,
z_thresh = z_thresh,
mean2use = mean2use,
sd2use = sd2use)
# mean2use = mean(test_data[,vars2use[1]] - test_data[,vars2use[2]])
# sd2use = sd(test_data[,vars2use[1]] - test_data[,vars2use[2]])
tmp_test = make_subtype(data2use = test_data,
z_thresh = z_thresh,
mean2use = mean2use,
sd2use = sd2use)
# compute model
mod2use = svm(x = tmp_train[,vars2use], y = tmp_train$z_ds_group)
pred_labels = predict(mod2use, tmp_test[,vars2use])
confmat = table(tmp_test$z_ds_group,pred_labels)
acc = (confmat[1,1]+confmat[2,2]+confmat[3,3])/length(pred_labels)
# plot confusion matrix
setHook("grid.newpage", function() pushViewport(viewport(x=1,y=1,width=0.9, height=0.9, name="vp", just=c("right","top"))), action="prepend")
pheatmap(confmat/rowSums(confmat)*100, display_numbers = TRUE, color = colorRampPalette(c('white','red'))(100), cluster_rows = FALSE, cluster_cols = FALSE, fontsize_number = 15,labels_row = c("RRB>SC","SC=RRB","SC>RRB"),labels_col = c("RRB>SC","SC=RRB","SC>RRB"),angle_col=90)
setHook("grid.newpage", NULL, "replace")
grid::grid.text("Actual Labels", y=-0.07, gp=gpar(fontsize=16))
grid::grid.text("Predicted Labels", x=-0.07, rot=90, gp=gpar(fontsize=16))
# show accuracy
true_accuracy = acc
true_accuracy
## [1] 0.958427
nperm = 10000
# make subtypes using z-scores computed from the mean and sd of the training set
train_data = Dverbal_Discovery
test_data = Dverbal_Replication
mean2use = mean(train_data[,vars2use[1]] - train_data[,vars2use[2]])
sd2use = sd(train_data[,vars2use[1]] - train_data[,vars2use[2]])
tmp_train = make_subtype(data2use = train_data,
z_thresh = z_thresh,
mean2use = mean2use,
sd2use = sd2use)
# mean2use = mean(test_data[,vars2use[1]] - test_data[,vars2use[2]])
# sd2use = sd(test_data[,vars2use[1]] - test_data[,vars2use[2]])
tmp_test = make_subtype(data2use = test_data,
z_thresh = z_thresh,
mean2use = mean2use,
sd2use = sd2use)
acc = vector(length = nperm)
for (iperm in 1:nperm){
# set seed for reproducibility
set.seed(iperm)
# compute model
permuted_labels = sample(tmp_train$z_ds_group)
mod2use = svm(x = tmp_train[,vars2use], y = permuted_labels)
pred_labels = predict(mod2use, tmp_test[,vars2use])
confmat = table(tmp_test$z_ds_group,pred_labels)
acc[iperm] = (confmat[1,1]+confmat[2,2]+confmat[3,3])/length(pred_labels)
} # for (iperm in 1:nperm)
df2plot = data.frame(Accuracy = acc)
p = ggplot(data = df2plot, aes(x = Accuracy)) + geom_histogram() + geom_vline(xintercept=true_accuracy)
p
# compute p-value
pval = sum(c(true_accuracy,acc)>=true_accuracy)/(nperm+1)
pval
## [1] 9.999e-05
maxScores = c(3,4)
# Discovery - make plot with all individuals shown as lines
df2use = Dverbal_Discovery[,c("subjectkey",adi_total_vars2use)]
df2use$subgrp = factor(tmp_train$z_ds_group)
df2use = data.frame(df2use)
df2use$subjectkey = factor(df2use$subjectkey)
df2use$SC = df2use$dbaes_atotal
df2use$RRB = df2use$dbaes_btotal
df4plot = melt(df2use,
id.vars = c("subjectkey","subgrp"),
measure.vars = c("SC","RRB"))
p = ggplot(data = df4plot, aes(x = variable,
y = value,
colour = subgrp,
group = subjectkey)) + facet_grid(. ~ subgrp)
p = p + geom_point(shape=1) + geom_line(alpha = 0.2) + ylim(0,1) + guides(color=FALSE)
p = p + ylab("Percent Severity") + xlab("ADI-R subscale")
p3_middle_top = p + guides(colour=FALSE)
ggsave(filename = file.path(plotpath, sprintf("final_NDAR_Disc_jitterplot_z%s.pdf",as.character(z_thresh))), plot = p3_middle_top)
p
maxScores = c(3,4)
# Replication - make plot with all individuals shown as lines
df2use = Dverbal_Replication[,c("subjectkey",adi_total_vars2use)]
df2use$subgrp = factor(tmp_test$z_ds_group)
df2use = data.frame(df2use)
df2use$subjectkey = factor(df2use$subjectkey)
df2use$SC = df2use$dbaes_atotal
df2use$RRB = df2use$dbaes_btotal
df4plot = melt(df2use,
id.vars = c("subjectkey","subgrp"),
measure.vars = c("SC","RRB"))
p = ggplot(data = df4plot, aes(x = variable,
y = value,
colour = subgrp,
group = subjectkey)) + facet_grid(. ~ subgrp)
p = p + geom_point(shape=1) + geom_line(alpha = 0.2) + ylim(0,1) + guides(color=FALSE)
p = p + ylab("Percent Severity") + xlab("ADI-R subscale")
p4_middle_bottom = p + guides(colour=FALSE)
ggsave(filename = file.path(plotpath, sprintf("final_NDAR_Rep_jitterplot_z%s.pdf",as.character(z_thresh))), plot = p4_middle_bottom)
p
# make SC and RRB percentages since EU-AIMS data is specified as percentages
Dverbal = read.csv(file.path(datapath,"tidy_verbal.csv"))
euaims_data = read.csv(file.path(datapath,"tidy_euaims.csv"))
mask1 = euaims_data$Diagnosis=="ASD"
mask2 = (is.na(euaims_data$A1_pct_severity) | is.na(euaims_data$A2_pct_severity) | is.na(euaims_data$A3_pct_severity) | is.na(euaims_data$B1_pct_severity) | is.na(euaims_data$B2_pct_severity) | is.na(euaims_data$B3_pct_severity) | is.na(euaims_data$B4_pct_severity))
euaims_data = subset(euaims_data, (mask1 & !mask2))
Dverbal[,vars2use[1]] = Dverbal[,vars2use[1]]
Dverbal[,vars2use[2]] = Dverbal[,vars2use[2]]
euaims_data[,vars2use[1]] = (euaims_data$A1_pct_severity +
euaims_data$A2_pct_severity +
euaims_data$A3_pct_severity)/3
euaims_data[,vars2use[2]] = (euaims_data$B1_pct_severity +
euaims_data$B2_pct_severity +
euaims_data$B3_pct_severity +
euaims_data$B4_pct_severity)/4
train_data = Dverbal
test_data = euaims_data
mean2use = mean(train_data[,vars2use[1]] - train_data[,vars2use[2]], na.rm=TRUE)
sd2use = sd(train_data[,vars2use[1]] - train_data[,vars2use[2]], na.rm=TRUE)
c(mean2use, sd2use)
## [1] -0.01045243 0.19482749
tmp_train = make_subtype(data2use = train_data,
z_thresh = z_thresh,
mean2use = mean2use,
sd2use = sd2use)
tmp_test = make_subtype(data2use = test_data,
z_thresh = z_thresh,
mean2use = mean2use,
sd2use = sd2use)
# compute model
mod2use = svm(x = tmp_train[,vars2use], y = tmp_train$z_ds_group)
pred_labels = predict(mod2use, tmp_test[,vars2use])
confmat = table(tmp_test$z_ds_group,pred_labels)
acc = (confmat[1,1]+confmat[2,2]+confmat[3,3])/length(pred_labels)
tmp_test$svm_pred_labels = pred_labels
# plot confusion matrix
setHook("grid.newpage", function() pushViewport(viewport(x=1,y=1,width=0.9, height=0.9, name="vp", just=c("right","top"))), action="prepend")
pheatmap(confmat/rowSums(confmat)*100, display_numbers = TRUE, color = colorRampPalette(c('white','red'))(100), cluster_rows = FALSE, cluster_cols = FALSE, fontsize_number = 15,labels_row = c("RRB>SC","SC=RRB","SC>RRB"),labels_col = c("RRB>SC","SC=RRB","SC>RRB"),angle_col=90)
setHook("grid.newpage", NULL, "replace")
grid::grid.text("Actual Labels", y=-0.07, gp=gpar(fontsize=16))
grid::grid.text("Predicted Labels", x=-0.07, rot=90, gp=gpar(fontsize=16))
# scatterplot
p1 = ggplot(data = tmp_train, aes(x = dbaes_atotal, y = dbaes_btotal, colour = factor(z_ds_group))) + geom_point() + xlab("Social-Communication") + ylab("Restricted Repetitive Behaviors") + ylim(0,0.8) + xlim(0,0.8) + ggtitle("NDAR ALL")
p1
p2 = ggplot(data = tmp_test, aes(x = dbaes_atotal, y = dbaes_btotal, colour = factor(z_ds_group))) + geom_point() + xlab("Social-Communication") + ylab("Restricted Repetitive Behaviors") + ylim(0,0.8) + xlim(0,0.8) + ggtitle("EU-AIMS with Groups from NDAR ALL")
p2
p3 = ggplot(data = tmp_test, aes(x = dbaes_atotal, y = dbaes_btotal, colour = factor(pred_labels))) + geom_point() + xlab("SC") + ylab("RRB") + ylim(0,0.8) + xlim(0,0.8) + ggtitle("EU-AIMS")
p5_bottom_right = p3 + guides(colour=FALSE)
ggsave(filename = file.path(plotpath, sprintf("final_EUAIMS_scatterplot_z%s.pdf",as.character(z_thresh))), plot = p5_bottom_right)
p3
# # write out EU-AIMS LEAP data with subgroups defined by NDAR All
write.csv(tmp_test, file.path(datapath, sprintf("tidy_euaims_NDAR_subtypes_diffscore_z%s.csv",as.character(z_thresh))))
Make final plot
p_final = p1_top_left + p3_middle_top + plot_spacer() + p2_bottom_left + p4_middle_bottom + p5_bottom_right + plot_layout(nrow=3, ncol=3, widths = c(4,4,4), heights = c(8,8,8))
ggsave(filename = file.path(plotpath, sprintf("final_NDAR_EUAIMS_subtypes_plot_z%s.pdf",as.character(z_thresh))), plot = p_final)
p_final
Integrate subtypes with rest of EU-AIMS LEAP data
euaims_data = read.csv(file.path(datapath,"tidy_euaims.csv"))
td_df = subset(euaims_data, euaims_data$Diagnosis=="TD",select=2:6)
td_df$A_pct_severity = as.numeric(NA)
td_df$B_pct_severity = as.numeric(NA)
td_df$subgrp = "TD"
cols2use = c("subid","age","Centre","Schedule","Diagnosis","dbaes_atotal","dbaes_btotal","svm_pred_labels")
tmp_asd = tmp_test[,cols2use]
colnames(tmp_asd)[6] = "A_pct_severity"
colnames(tmp_asd)[7] = "B_pct_severity"
colnames(tmp_asd)[8] = "subgrp"
asd_df = tmp_asd
all_data = rbind(td_df,asd_df)
fname = "/Users/mlombardo/Dropbox/euaims/data/rsfmri_preproc/euaims_preproc.xlsx"
pp_data = read_excel(fname)
mask = pp_data$notes=="ok"
pp_data = subset(pp_data, mask)
asd_df = merge(pp_data[,c("subid","sex")],asd_df, by = "subid")
td_df = merge(pp_data[,c("subid","sex")],td_df, by = "subid")
all_data = rbind(td_df,asd_df)
data2write = merge(pp_data, all_data, by = "subid")
data2write$age = data2write$age.x
data2write$sex = data2write$sex.y
print(table(data2write$Schedule, data2write$subgrp))
##
## RRB_over_SC SC_equal_RRB SC_over_RRB TD
## A 6 39 35 78
## B 1 51 33 83
## C 4 33 26 59
## D 0 13 25 23
print(table(data2write$Centre, data2write$subgrp))
##
## RRB_over_SC SC_equal_RRB SC_over_RRB TD
## CAMBRIDGE 4 30 8 29
## KINGS_COLLEGE 6 51 49 78
## MANNHEIM 0 0 0 34
## NIJMEGEN 1 43 49 64
## UTRECHT 0 12 13 38
print(table(data2write$sex, data2write$subgrp))
##
## RRB_over_SC SC_equal_RRB SC_over_RRB TD
## Female 5 35 31 88
## Male 6 101 88 155
#DSM-5 - find best split that balances participants across sites
a = findBestSplit(asd_df, seed_range = c(172342)) #,300001:500000))
## [1] 172342
best_seeds = a$seed[!is.na(a$discrepancy) & a$discrepancy==min(a$discrepancy, na.rm = TRUE)]
print(best_seeds)
## [1] 172342
# Split datasets -------------------------------------------------------------
rngSeed = best_seeds[1]
# split Schedule A dataset
dset2use = subset(asd_df, asd_df$Schedule=="A")
tmp_d = SplitDatasetsBySex(dset2use, rngSeed = rngSeed)
A_Discovery = tmp_d[[2]]
A_Replication = tmp_d[[1]]
# split Schedule B dataset
dset2use = subset(asd_df, asd_df$Schedule=="B")
tmp_d = SplitDatasetsBySex(dset2use, rngSeed = rngSeed)
B_Discovery = tmp_d[[2]]
B_Replication = tmp_d[[1]]
# split Schedule C dataset
dset2use = subset(asd_df, asd_df$Schedule=="C")
tmp_d = SplitDatasetsBySex(dset2use, rngSeed = rngSeed)
C_Discovery = tmp_d[[2]]
C_Replication = tmp_d[[1]]
# split Schedule D dataset
dset2use = subset(asd_df, asd_df$Schedule=="D")
tmp_d = SplitDatasetsBySex(dset2use, rngSeed = rngSeed)
D_Discovery = tmp_d[[2]]
D_Replication = tmp_d[[1]]
df_Disc = rbind(A_Discovery, B_Discovery, C_Discovery, D_Discovery)
df_Rep = rbind(A_Replication, B_Replication, C_Replication, D_Replication)
a = table(df_Disc$Schedule, df_Disc$Centre); print(a)
##
## CAMBRIDGE KINGS_COLLEGE NIJMEGEN UTRECHT
## A 6 17 10 7
## B 9 16 14 3
## C 6 9 14 2
## D 0 10 10 0
b = table(df_Rep$Schedule, df_Rep$Centre); print(b)
##
## CAMBRIDGE KINGS_COLLEGE NIJMEGEN UTRECHT
## A 7 18 10 5
## B 8 17 13 5
## C 6 10 13 3
## D 0 9 9 0
print(a-b)
##
## CAMBRIDGE KINGS_COLLEGE NIJMEGEN UTRECHT
## A -1 -1 0 2
## B 1 -1 1 -2
## C 0 -1 1 -1
## D 0 1 1 0
print(sum(rowSums(abs(a-b))))
## [1] 14
data2write$dataset = NA
mask = is.element(data2write$subid,df_Disc$subid)
data2write[mask,"dataset"] = "Discovery"
mask = is.element(data2write$subid,df_Rep$subid)
data2write[mask,"dataset"] = "Replication"
asd_Disc = subset(data2write, data2write$dataset=="Discovery" & data2write$Diagnosis=="ASD")
asd_Rep = subset(data2write, data2write$dataset=="Replication" & data2write$Diagnosis=="ASD")
# find which seed gives best TD age-match -------------------------------------
seeds = 1:1000
pvals = data.frame(matrix(nrow = length(seeds), ncol = 2))
for (i in 1:length(seeds)) {
res = findTDAgeMatch(data2write, seed_range = c(seeds[i],seeds[i]))
td_Disc_matched = res[[2]]
td_Rep_matched = res[[1]]
tres = t.test(td_Disc_matched$age, asd_Disc$age)
pvals[i,1] = tres$p.value
tres = t.test(td_Rep_matched$age, asd_Rep$age)
pvals[i,2] = tres$p.value
#print(i)
}
a = sort.int(pvals[,1], decreasing = TRUE, index.return = TRUE)
b=pvals[a$ix,]
seed2use = 929
res = findTDAgeMatch(data2write, seed_range = c(seed2use,seed2use))
td_Disc_matched = res[[2]]
td_Rep_matched = res[[1]]
mask = is.element(data2write$subid, td_Disc_matched$subid)
data2write$dataset[mask] = "Discovery"
mask = is.element(data2write$subid, td_Rep_matched$subid)
data2write$dataset[mask] = "Replication"
print(table(data2write$dataset, data2write$subgrp))
##
## RRB_over_SC SC_equal_RRB SC_over_RRB TD
## Discovery 8 66 59 121
## Replication 3 70 60 122
fname2save = here(sprintf("asd_subgrp_data_rsfmri_ALL_DSM5_diffzscoreGrps_z%s.csv",as.character(z_thresh)))
write.csv(data2write,file = fname2save)
Plot mean FD and hypothesis test
df2use = subset(data2write, data2write$subgrp!="RRB_over_SC")
p = ggplot(data = df2use, aes(x = subgrp, y =meanFD, colour = subgrp)) + facet_wrap(. ~ dataset)
p = p + geom_jitter() + geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA) + guides(colour = FALSE)
p = p + xlab("Group") + ylab("Mean FD")
p
# Hypothesis test
df2use = data2write
y_var = "meanFD"
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df2use, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 501 116.11403 <.0001
## subgrp 3 501 1.94055 0.1221
# Discovery
df4mod = subset(df2use,df2use$dataset=="Discovery" & df2use$subgrp!="RRBoverSC")
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df4mod, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 246 145.64335 <.0001
## subgrp 3 246 2.73468 0.0443
# Replication
df4mod = subset(df2use,df2use$dataset=="Replication" & df2use$subgrp!="RRBoverSC")
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df4mod, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 247 67.03279 <.0001
## subgrp 3 247 0.21626 0.8851
Plot ADOS scores in the subtypes and also run a hypothesis test
df2use = subset(data2write, data2write$subgrp!="TD")
df2use$ados_2_SA_CSS[df2use$ados_2_SA_CSS==999] = NA
df2use$ados_2_RRB_CSS[df2use$ados_2_RRB_CSS==999] = NA
df4plot = melt(df2use,
id.vars = c("subid","subgrp"),
measure.vars = c("ados_2_SA_CSS","ados_2_RRB_CSS"))
df4plot$variable = as.character(df4plot$variable)
df4plot$variable[df4plot$variable=="ados_2_SA_CSS"] = "SA"
df4plot$variable[df4plot$variable=="ados_2_RRB_CSS"] = "RRB"
df4plot$variable = factor(df4plot$variable)
p = ggplot(data = df4plot, aes(x = variable,
y = value,
colour = subgrp,
group = subid)) + facet_grid(. ~ subgrp)
p = p + geom_point(shape=1) + geom_line(alpha = 0.2) + guides(color=FALSE)
p = p + ylab("ADOS CSS") + xlab("ADOS Subscale")
ggsave(filename = file.path(plotpath, sprintf("final_EUAIMS_ADOS_jitterplot_z%s.pdf",as.character(z_thresh))))
p
# hypothsis test on ADOS SA CSS
y_var = "ados_2_SA_CSS"
df2use = subset(data2write, !is.element(data2write$subgrp,c("TD")))
# Hypothesis test
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df2use, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 260 9.693338 0.0021
## subgrp 2 260 0.587368 0.5565
# hypothsis test on ADOS RRB CSS
y_var = "ados_2_RRB_CSS"
df2use = subset(data2write, !is.element(data2write$subgrp,c("TD")))
# Hypothesis test
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df2use, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 260 9.917288 0.0018
## subgrp 2 260 0.571781 0.5652
Descriptive stats
df2use = data2write[,c("subid","Diagnosis","dataset","Centre","meanFD","age","sex","subgrp","viq_all","piq_all","fsiq4_all","ados_2_SA_CSS","ados_2_RRB_CSS","SRS_tscore","SRS_tscore_self")]
mask = df2use==999
df2use[mask] = NA
df2use$age = df2use$age/365
cols2use = c("dataset","subgrp","age","meanFD","viq_all","piq_all","fsiq4_all","ados_2_SA_CSS","ados_2_RRB_CSS","SRS_tscore","SRS_tscore_self")
res=describeBy(x = df2use[,cols2use], group = c("subgrp"))
res
##
## Descriptive statistics by group
## group: RRB_over_SC
## vars n mean sd median trimmed mad min max range
## dataset* 1 11 NaN NA NA NaN NA Inf -Inf -Inf
## subgrp* 2 11 NaN NA NA NaN NA Inf -Inf -Inf
## age 3 11 16.72 6.37 19.56 16.90 6.41 7.89 23.88 16.00
## meanFD 4 11 0.26 0.31 0.18 0.18 0.10 0.06 1.14 1.09
## viq_all 5 11 105.80 17.44 106.00 104.76 14.83 78.00 143.00 65.00
## piq_all 6 11 106.27 17.60 100.00 103.78 16.31 87.00 148.00 61.00
## fsiq4_all 7 11 106.27 17.57 104.00 104.56 11.86 80.00 148.00 68.00
## ados_2_SA_CSS 8 11 4.00 2.68 3.00 3.78 2.97 1.00 9.00 8.00
## ados_2_RRB_CSS 9 11 4.00 3.61 1.00 3.78 0.00 1.00 9.00 8.00
## SRS_tscore 10 7 70.43 10.74 72.00 70.43 8.90 58.00 90.00 32.00
## SRS_tscore_self 11 5 60.60 7.13 63.00 60.60 5.93 49.00 67.00 18.00
## skew kurtosis se
## dataset* NA NA NA
## subgrp* NA NA NA
## age -0.17 -1.96 1.92
## meanFD 2.13 3.39 0.09
## viq_all 0.52 -0.35 5.26
## piq_all 1.02 0.18 5.31
## fsiq4_all 0.88 0.42 5.30
## ados_2_SA_CSS 0.51 -1.27 0.81
## ados_2_RRB_CSS 0.35 -1.88 1.09
## SRS_tscore 0.51 -1.03 4.06
## SRS_tscore_self -0.66 -1.43 3.19
## ------------------------------------------------------------
## group: SC_equal_RRB
## vars n mean sd median trimmed mad min max range
## dataset* 1 136 NaN NA NA NaN NA Inf -Inf -Inf
## subgrp* 2 136 NaN NA NA NaN NA Inf -Inf -Inf
## age 3 136 16.27 5.73 15.49 15.94 5.60 7.08 30.28 23.20
## meanFD 4 136 0.28 0.40 0.19 0.22 0.12 0.04 3.95 3.91
## viq_all 5 134 99.86 17.37 102.00 100.28 16.31 61.00 136.00 75.00
## piq_all 6 134 102.04 20.07 103.49 102.80 19.97 52.00 150.00 98.00
## fsiq4_all 7 135 101.18 17.75 104.00 101.91 19.27 60.00 143.00 83.00
## ados_2_SA_CSS 8 133 5.98 2.48 6.00 6.07 2.97 1.00 10.00 9.00
## ados_2_RRB_CSS 9 133 4.98 2.67 6.00 4.98 1.48 1.00 10.00 9.00
## SRS_tscore 10 121 69.19 11.97 69.00 69.26 13.34 43.00 90.00 47.00
## SRS_tscore_self 11 65 61.69 9.79 62.00 61.43 10.38 43.00 89.00 46.00
## skew kurtosis se
## dataset* NA NA NA
## subgrp* NA NA NA
## age 0.49 -0.48 0.49
## meanFD 6.46 52.97 0.03
## viq_all -0.22 -0.70 1.50
## piq_all -0.31 -0.29 1.73
## fsiq4_all -0.31 -0.55 1.53
## ados_2_SA_CSS -0.29 -0.88 0.22
## ados_2_RRB_CSS -0.35 -0.94 0.23
## SRS_tscore -0.03 -0.91 1.09
## SRS_tscore_self 0.26 -0.23 1.21
## ------------------------------------------------------------
## group: SC_over_RRB
## vars n mean sd median trimmed mad min max range
## dataset* 1 119 NaN NA NA NaN NA Inf -Inf -Inf
## subgrp* 2 119 NaN NA NA NaN NA Inf -Inf -Inf
## age 3 119 16.47 5.15 16.14 16.32 6.02 7.48 29.40 21.92
## meanFD 4 119 0.22 0.21 0.15 0.18 0.10 0.03 1.31 1.27
## viq_all 5 115 96.24 19.51 98.00 96.60 20.76 50.91 142.00 91.09
## piq_all 6 117 98.16 21.33 103.00 99.49 20.76 44.03 138.00 93.97
## fsiq4_all 7 117 97.63 19.30 102.00 98.30 19.92 59.00 139.00 80.00
## ados_2_SA_CSS 8 114 6.22 2.79 7.00 6.34 2.97 1.00 10.00 9.00
## ados_2_RRB_CSS 9 114 4.59 2.75 5.00 4.48 2.97 1.00 10.00 9.00
## SRS_tscore 10 104 74.05 11.66 76.00 74.81 11.86 44.00 95.00 51.00
## SRS_tscore_self 11 53 63.11 11.03 61.00 62.30 8.90 40.00 94.00 54.00
## skew kurtosis se
## dataset* NA NA NA
## subgrp* NA NA NA
## age 0.29 -0.63 0.47
## meanFD 2.78 8.87 0.02
## viq_all -0.16 -0.45 1.82
## piq_all -0.52 -0.58 1.97
## fsiq4_all -0.33 -0.89 1.78
## ados_2_SA_CSS -0.30 -1.12 0.26
## ados_2_RRB_CSS -0.14 -1.17 0.26
## SRS_tscore -0.48 -0.46 1.14
## SRS_tscore_self 0.68 0.66 1.51
## ------------------------------------------------------------
## group: TD
## vars n mean sd median trimmed mad min max range
## dataset* 1 243 NaN NA NA NaN NA Inf -Inf -Inf
## subgrp* 2 243 NaN NA NA NaN NA Inf -Inf -Inf
## age 3 243 16.84 5.66 16.58 16.63 6.57 6.89 29.84 22.95
## meanFD 4 243 0.20 0.34 0.13 0.15 0.07 0.03 4.60 4.58
## viq_all 5 241 104.27 18.62 106.00 105.33 14.83 45.00 160.00 115.00
## piq_all 6 241 105.36 18.92 107.00 107.07 16.31 49.00 147.00 98.00
## fsiq4_all 7 241 105.33 17.70 108.18 107.16 13.61 50.00 142.00 92.00
## ados_2_SA_CSS 8 0 NaN NA NA NaN NA Inf -Inf -Inf
## ados_2_RRB_CSS 9 0 NaN NA NA NaN NA Inf -Inf -Inf
## SRS_tscore 10 133 47.54 9.34 44.00 45.82 4.45 37.00 90.00 53.00
## SRS_tscore_self 11 132 47.50 5.90 46.00 46.72 5.19 39.00 69.00 30.00
## skew kurtosis se
## dataset* NA NA NA
## subgrp* NA NA NA
## age 0.28 -0.74 0.36
## meanFD 9.47 111.25 0.02
## viq_all -0.55 0.92 1.20
## piq_all -0.80 0.53 1.22
## fsiq4_all -0.96 1.13 1.14
## ados_2_SA_CSS NA NA NA
## ados_2_RRB_CSS NA NA NA
## SRS_tscore 1.85 3.60 0.81
## SRS_tscore_self 1.18 1.39 0.51
Model age differences
y_var = "age"
# df4mod = subset(df2use,df2use$subgrp!="RRBoverSC")
df4mod = df2use
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df4mod, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 501 788.6982 <.0001
## subgrp 3 501 0.5741 0.6323
# Discovery
df4mod = subset(df2use,df2use$dataset=="Discovery" & df2use$subgrp!="RRBoverSC")
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df4mod, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 246 600.5126 <.0001
## subgrp 3 246 0.3079 0.8197
# Replication
df4mod = subset(df2use,df2use$dataset=="Replication" & df2use$subgrp!="RRBoverSC")
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df4mod, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 247 484.5865 <.0001
## subgrp 3 247 0.5409 0.6547
Model VIQ differences
y_var = "viq_all"
# df4mod = subset(df2use,df2use$subgrp!="RRBoverSC")
df4mod = df2use
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df4mod, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 493 2033.3530 <.0001
## subgrp 3 493 4.0341 0.0075
# Discovery
df4mod = subset(df2use,df2use$dataset=="Discovery" & df2use$subgrp!="RRBoverSC")
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df4mod, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 242 2242.2320 <.0001
## subgrp 3 242 2.4312 0.0658
# Replication
df4mod = subset(df2use,df2use$dataset=="Replication" & df2use$subgrp!="RRBoverSC")
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df4mod, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 243 1575.8466 <.0001
## subgrp 3 243 2.8377 0.0387
Model PIQ differences
y_var = "piq_all"
# df4mod = subset(df2use,df2use$subgrp!="RRBoverSC")
df4mod = df2use
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df4mod, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 495 1495.1212 <.0001
## subgrp 3 495 1.7613 0.1536
# Discovery
df4mod = subset(df2use,df2use$dataset=="Discovery" & df2use$subgrp!="RRBoverSC")
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df4mod, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 242 1868.657 <.0001
## subgrp 3 242 1.567 0.198
# Replication
df4mod = subset(df2use,df2use$dataset=="Replication" & df2use$subgrp!="RRBoverSC")
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df4mod, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 245 1228.873 <.0001
## subgrp 3 245 1.279 0.2821
Model FIQ differences
y_var = "fsiq4_all"
# df4mod = subset(df2use,df2use$subgrp!="RRBoverSC")
df4mod = df2use
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df4mod, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 496 2110.521 <.0001
## subgrp 3 496 3.060 0.0279
# Discovery
df4mod = subset(df2use,df2use$dataset=="Discovery" & df2use$subgrp!="RRBoverSC")
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df4mod, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 244 2782.0395 <.0001
## subgrp 3 244 2.3548 0.0726
# Replication
df4mod = subset(df2use,df2use$dataset=="Replication" & df2use$subgrp!="RRBoverSC")
# construct linear model
# mixed-effect model: site as random factor, all other covariates as fixed factors
fx_form = as.formula(sprintf("%s ~ %s",y_var,"subgrp"))
rx_form = as.formula(sprintf("~ 1|%s","Centre"))
mod2use = eval(substitute(lme(fixed = fx_form, random = rx_form, data = df4mod, na.action = na.omit)))
# run ANOVA
res = anova(mod2use)
res
## numDF denDF F-value p-value
## (Intercept) 1 244 1535.0793 <.0001
## subgrp 3 244 1.8337 0.1416